Exact linear hydrodynamics from the Boltzmann equation 
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Exact (to all orders in Knudsen number) equations of linear hydrodynamics are derived from the Boltzmann 
kinetic equation with the Bhatnagar-Gross-Krook collision integral. The exact hydrodynamic equations are cast 
in a form which allows us to immediately prove their hyperbolicity, stability, and existence of an //-theorem. 
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Hydrodynamics assumes that a state of a fluid is solely de- 
scribed by five fields: density, momentum and temperature. 
Derivation of the Navier-Stokes-Fourier (NSF) hydrodynamic 
equations from the Boltzmann kinetic equation as the first- 
order approximation in the Knudsen number (ratio between 
mean free path and a flow scale) by Enskog and Chapman 
is a textbook example of a success of statistical physics [1]. 
Recent renewed interest to the problems beyond the standard 
hydrodynamics is due, in particular, to flow simulation and ex- 
periments at a micro- and nano-scale [2—4]. However, almost 
a century of effort to extend the hydrodynamic description be- 
yond the NSF approximation failed even in the case of small 
deviations around the equilibrium. In order to appreciate the 
problem, let us remind that, in the NSF approximations, the 
decay rate of the hydrodynamic modes is quadratic in the 
wave vector, Rc(oj) ~ — fc 2 , and is unbounded. On the other 
hand, Boltzmann's collision term features equilibration with 
finite characteristic rates. This "finite collision frequency" is 
obviously incompatible with the arbitrary decay rates in the 
NSF approximation: intuitively, hydrodynamic modes at large 
k cannot relax faster than the collision frequency. Now, the 
classical method of Enskog and Chapman extends the hydro- 
dynamics beyond the NSF in such a way that the decay rate of 
the next order approximations (Burnett and super-Burnett) are 
polynomials of higher order in k. In such an extension, relax- 
ation rate may become completely unphysical (amplification 
instead of attenuation), as first shown by Bobylev [5] for a 
particular case of Maxwell molecules. This indicates inability 
of the Chapman-Enskog method to tackle the above problem, 
and non-perturbative approaches are sought. The problem of 
exact hydrodynamics has been studied in depth recently for 
toy (finite-dimensional) models - moment systems of Grad - in 
[6-8], and many remarkable results were obtained. In partic- 
ular, in [6, 7] it was shown that the exact hydrodynamic equa- 
tions are hyperbolic and stable for all wave numbers. How- 
ever, for "true" kinetic equations such questions remain open. 

In this Letter we derive exact hydrodynamic equations from 
the linearized Boltzmann equation with the Bhatnagar-Gross- 
Krook (BGK) collision term. This kinetic equation remains 
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FIG. 1: Exact hydrodynamic modes w of the Boltzmann-BGK kinetic 
equation as a function of wave number k (two complex-conjugated acous- 
tic modes oj a c , twice degenerated shear mode aj s h C ar and thermal diffusion 
mode t^dift )- The non-positive decay rates Rc(o;) attain the limit of collision 
frequency (— 1) as fc oo. 



popular in applications [9], and features a single relaxation 
rate. The result for the hydrodynamic modes is demonstrated 
in Fig. 1. It is clear from Fig. 1 that the relaxation of none of 
the hydrodynamic modes is faster than uj = — 1 which is the 
collision frequency in the units adopted in this paper. Thus, 
the result for the exact hydrodynamics indeed corresponds to 
the above intuitive picture. Below, we apply the method of in- 
variant manifold [10] to derive the hydrodynamic equations. 
The non-perturbative derivation is made possible with an op- 
timal combination of analytical and numerical approaches to 
solve the invariance equation. 

Point of departure is the linearized Boltzmann-BGK equa- 
tion for the deviation A/ = / — / GM of the distribution func- 
tion / from a global Maxwellian / GM (c 2 ) = 7r~ 3 / 2 e~ c . In 
the reciprocal space, it reads, 



d t Af = -ik-cAf-Sf; Sf = f-f 



LM 



(1) 
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with the wave vector k = k e y defining e y , k = |k|, peculiar 
velocity c and time t. All quantities are considered dimen- 
sionless, i.e., reduced with the un its of the relaxation time t, 
the thermal velocity y/2ksT '/m and mass m of the particle. 
In (1), the linearized local Maxwellian is / LM = f GM (l+ip ) 
where 1 + ip Q = + 2c ■ (c) f + |(c 2 - |)(c 2 - §),. Av- 
erages are defined for arbitrary u> via (u>) f = J d 3 c, and 
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we introduce pertinent quantities which characterize deviation 
from the global equilibrium: n = (l)f — 1 (density perturba- 
tion), u = (c) f (velocity perturbation) and T = | (c 2 — |)/ 
(temperature perturbation). Since the scalar product between 
k and c appears in (1), the distribution function offers symme- 
try with respect to the ey-axis, which is not uniaxial in case u 
is not collinear with e B . We denote the two components of 
the mean velocity as u" = u • e B and u x = u ■ e ± , where 
the unit vector e ± belongs to the intersection of the plane per- 
pendicular to k and the plane spanned by k and u, so that 
u = u"e|| + w i e i . Equations of change for moments (a;) are 
obtained by integration of the weighted (1) over d 3 c as 

9t(w) A/ = -ik-(cw) A/ -(w) (5/ . (2) 

In order to calculate such averages, we can switch to spher- 
ical coordinates. For each (at present arbitrary) wave vec- 
tor, we choose the coordinate system in such a way that its 
z-direction aligns with e B . We can then express c in terms 
of its norm c, a vertical variable z and plane vector e^, (az- 
imuthal angle ■ e ± = cos 0) for the present purpose as 
c/c = Vl — z 2 + ze y . We could have equally chosen a 
fixed coordinate system in the plane orthogonal to k, and two 
fields instead of m 1 plus an angle, viz. u^+tUy = e^u^. Due 
to isotropy, u 1 alone fully represents the twice degenerated 
(shear) dynamics. In order to simplify notation and compute 
the dynamics of all five fields we introduce a four-dimensional 
vector of hydrodynamic fields, x = [xi,X2,x^,Xi) = 
(x !l , X4) with x" = (n, u" , T) and X4 = m 1 . Then ipo takes a 
simple form, (po = X° ■ x. The vector X° can immediately be 
read off, we have X°(c, z) = (1, 2c B , c 2 — §, 2c^), where we 
introduced, for later use, the abbreviations 

c,| = c • e,| , c,p = c-e ± , c ± = — , (3) 

e_L • e 

such that ik • c = ifccy, c = c ± e<p + Cyey with c y = cz and 
c x = cy/1 — z 2 , contrasted by (and e^), do not depend 
on the azimuthal angle. Similarly, we introduce yet unknown 
fields <5X(c, k) which characterize the nonequilibrium part of 
the distribution function, Sip = Sf / f GM in terms of the hy- 
drodynamic fields x themselves, 

5tp = SX • x = 6X1 n + 5X 2 u 11 + SX 3 T + SX^, (4) 

where 6X4 factorizes as 6X4(0, z, 0) = 26Yi(c, z) e# • e_i_. 
This "eigen"-closure (4) which formally and very generally 
addresses the fact, that we wish to not include other than hy- 
drodynamic variables implies a closure between moments of 
the distribution function, to be worked out in detail below. It 
assumes the existence of an invariant manifold, and the hy- 
drodynamic fields as slow variables which leave the higher 
moments "slaved". In order to ensure these contributions to 
not interfere with the local Maxwellian, one has the freedom 
to require (X°)sf = without producing any limitation, i.e., 
by keeping x to be defined through the local Maxwellian part 
of the distribution function. Using the above form for Sf in 
(1), and using the canonical abbreviation AX = X° + <5X, 
yields 

AX ■ <9 t x = -ik c„ AX • x - SX ■ x, (5) 
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-k 2 B ikA -k 2 C ikD 
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real, © imag,© real,© imag,Q 
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(7»MTi) (t 11 ^) (7"c||5X 3 ) <(c 2 - l)c ± SY 4 ) 

ikX -k 2 Z ikY -k 2 U 
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TABLE I: Symmetry adapted components of (nonequilibrium) stress 
tensor cr and heat flux q, introduced in (7a) and (7b), respectively. 
Row 2: Microscopic expression of these components (averaging with 
the global Maxwellian). Short-hand notation used: A" = c 2 — ^- and 
7" = (c 2 — §)c||. Row 3: Expression of the components in terms 
of (as we show, real-valued) functions A-Z (see text). Row 4: Val- 
ues of functions A-Z at k — 0. These values recover hydrodynamic 
equations up to Burnett approximation. Row 5: Parity with respect 
to z - symmetric (©) or antisymmetric (©) - of the part of the cor- 
responding 8X entering the averaging in row 2, and whether this 
part is imaginary or real-valued (see Fig. 2). Row 3 is an immediate 
consequence of row 5. 

which is a nonlinear integral equation for the unknown fields 
SX, because <9 t x has to be replaced by the right hand side 
of (2), for a suitably chosen vector u> fulfilling (w)a/ = x. 
Here, u> is similar with X° and differs from X° mainly be- 
cause of conventions for prefactors in the temperature and ve- 
locity definitions, u; = (1, C|| , |(c 2 — |), Ctj,). Within the same 
eigen-closure, Eq. (2) is linear in x and hence written as 

<9 t x = M • x. (6) 

The matrix M solely depends on the non-hydrodynamic 
fields, the heat flux q = (c(c 2 — §))/ and the stress tensor 
<7 = ( cc 1 ) f, where s 1 denotes the symmetric traceless part 
of a tensor s [6, 11], s 1 = ±(s + s T ) - |tr(s)I. Using (4), 
the stress tensor and heat flux uniquely decompose as follows 

3 

a = ^ 2 e n e ii +^2 6,1^, (7a) 
q = q" e„ + q ± e ± , (7b) 

with the moments cr" = (<r\, a\, 03) ■ x" and = 04X4, and 
similarly for q (see Row 2 of Tab. I). The prefactors arise from 
the identities 'eyey' : e n e n — | and "e n e,| ' : e^e ± = |. The 
appearance of 8Y4 rather than SX4 in the expression for the 
orthogonal moment (in Tab. I) reflects the fact that we have al- 
ready integrated out the angular variable, J Q 27r e^e^ • e ± d<fi — 
ire ± . We note in passing that, while the stress tensor has, in 
general, three different eigenvalues, in the present symmetry 
adapted coordinate system it exhibits a vanishing first normal 
stress difference. Since the integral kernels of all moments in 
(7) do not depend on the azimuthal angle, these are actually 
two-dimensional integrals over c G [0, 00] and z e [—1,1], 
weighted by 27rc 2 / GM (c 2 )<5X Al . 
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Stress tensor and heat flux can yet be written in an alterna- 
tive form which is defined by Row 3 of Tab. I. As we will 
see later on, due to basic symmetry considerations, the hereby 
introduced functions A-Z are real-valued. We postpone the 
related proof, and proceed by using these functions A-Z to 
split M into parts as M = Re(M) - i Im(M) with 



Re(M) = fc 2 



Im(M) 
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FIG. 2: (Color online) Sample distribution function /(c, k) at k = 1, fully 
characterized by the four quantities SX± t 2,3(c, z) and 6Yi(c, z). Shown here 
are both their real (left) and imaginary parts (right column). In order to im- 
prove contrast, we actually plot In |1 + f GM SX l _ l \ multiplied by the sign of 
SXfj,. Same color code for all plots, ranging from —0.2 (red) to +0.2 (blue). 



Note that the checkerboard structure of the matrix M (8) 
is particularly useful for studying properties of the hydrody- 
namic equations (6), such as hyperbolicity and stability (see 



[6, 7] and below), once the functions A-Z are explicitly eval- 
uated. For that, we still require <5X. Combining (5) and (6), 
and requiring that the result holds for any x (invariance con- 
dition), we obtain a closed, singular integral equation (invari- 
ance equation) for complex-valued <5X, 



SX. = X° • (M + [ik C]] + I}!)' 1 - X°. 



(9) 



Notice that SX vanishes for k = 0, and that (9) is supple- 
mented with the basic constraint (X )^ = 0, or equally, 
vanishing Lagrange multipliers (matrix) (X°#X), which, 
however, is automatically dealt with if we only evaluate 
anisotropic (irreducible) moments with Sf, such as those 
listed in Tab. I. The implicit equation (9) is identical with 
the eigen-closure (4), and is our main and practically useful 
result, with M from (8), Cy, X°, and A-Z defined in and just 
before (3) and Tab. I, respectively. 

We iteratively calculate (i) SX directly from (9) for each k 
in terms of M, (ii) subsequently calculate moments from <5X 
by numerical integration. Importantly, the fix point of the it- 
eration (i)-(ii)-(i)-.. is unique for each k, i.e., does not depend 
on the initial values for moments A-Z, as long as we choose 
real-valued ones which are consistent with (9), as we prove 
in the next paragraph. In addition, two other computational 
strategies were implemented: First, we used continuation of 
functions A-Z from their values at k = to solve (9) with an 
incremental increase of k, where the solution at k was used as 
the initial guess for k + dk. Second, we used also a continu- 
ation "backwards" in which the solution at some k (obtained 
by convergent iterations with a random initial condition) was 
used as the initial guess for a solution at k — dk. Both these 
strategies returned the same values of functions A-Z as com- 
puted by iterations from arbitrary initial condition. The solu- 
tion <5X allows to calculate the whole distribution function / 
via (4) as illustrated by Fig. 2. For resulting moments for a 
wide range of fc-values see Fig. 3. 

Finally, we need to clarify the origin of row 5 in Tab. I 
(which is directly illustrated by Fig. 2) and its implication on 
the structure of M (8) whose entries are - a priori - complex- 
valued functions to be calculated with complex-valued <5X. 
We wish to make use of the fact that all integrals over z van- 
ish for odd integrands. To this end we introduce abbrevia- 
tions © (0) for a real-valued quantity which is even (odd) 
with respect to the transformation z — > —z. One notices 
X° = (ffi, e, ffi, ffi), and we recall that A-Z are integrals 
over either even or odd functions in z, times a component of 
SX (see Tab. I). Let us prove the consistency of the spec- 
ified symmetry of M and the invariance condition: Start by 
assuming A-Z to be real-valued functions. Then M^ IV = © 
if /i + v is even, and A/ M „ = iQ otherwise. This implies 
8X X = © + iQ, SX 2 = © + i®, 5X 3 = © + ie, and 
SX4 = © + iQ, i.e., different symmetry properties for real 
and imaginary parts. With these "symmetry" expressions for 
X°, (5X, and M at hand, we can insert into the right hand side 
of the equation, SX = (X° + SX) ■ (M + i © I), which is 
identical with the invariance equation (9). There are only two 
cases to consider, because M has a checkerboard structure, 
i.e., only two types of columns: Columns ji = 1 and fj, = 3: 



SX U 



iQ because M\^n = 0; Columns fx G {2, 4}: 
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SX^ = © + iO if M M .i_3 = (which is the case for column 
4) and + i© if M M 4 = (which is the case for column 2). 
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FIG. 3: Moments A—Z vs. wave number k obtained with the solution of (9). 



d t x' = M' • x\ with M' = T ■ M ■ T" 1 is manifestly hyper- 
bolic and stable; Im(M') is symmetric, Rc(M') is symmetric 
and non-positive semi-definite. The corresponding transfor- 
mation matrix T can be easily read off the results obtained in 
[6, 7] for Grad's systems since the structure of the matrix M 
(8) is identical to the one studied in [6, 7]. We have explic- 
itly verified that matrix T (equations (21)-(23) in Ref. [7] and 
(13) in Ref. [6]) with the functions A—Z derived herein is real- 
valued and thus render the transformed hydrodynamic equa- 
tions manifestly hyperbolic and stable. We note that this result 
- hyperbolicity of exact hydrodynamic equations - strongly 
supports a recent suggestion by Bobylev to consider a hyper- 
bolic regularization of the Burnett approximation [12]. Sim- 
ilarly, using the hyperbolicity, an iJ-theorem is elementary 
proven as in [6, 12]. Finally, using the accurate data for func- 
tions A—Z, we can write analytic approximations for the hy- 
drodynamic equations (6) in such a way that hyperbolicity and 
stability is not destroyed in such an approximation (see [7]). 



We have thus shown that both sides of the invariance equa- 
tion (9) have equal symmetry properties, and that <5X with the 
specified symmetries is consistent with real-valued moments 
A-Z. The proof implies, that any iteratively obtained solu- 
tion, if it exists, starting with arbitrary real-valued moments 
A-Z in (9) to evaluate <5X must converge to real-valued so- 
lution A-Z. Since the solution is smoothly varying with k, 
and since A-Z at k = are known and are real-valued, the 
moments must be real-valued over the whole fc-space. 

With the result for the functions A-Z at hand, the extended 
hydrodynamic equations are closed. Let us briefly discuss 
the pertinent properties of this system. First, the general- 
ized transport coefficients are given by the nontrivial eigen- 
values of — fc~ 2 Rc(M): A2 = —A (elongation viscosity), 
A3 = — I Y (thermal diffusivity), and A4 = — D (shear vis- 
cosity). All these generalized transport coefficients are non- 
negative (see Fig. 3). Second, computing the eigen-values 
of matrix M we obtain the dispersion relation oj(k) of the 
corresponding hydrodynamic modes already presented in Fig. 
1. Third, a suitable transform of the hydrodynamic fields, 
x' = T • x, where T is a real-valued matrix, can be es- 
tablished such that the transformed hydrodynamic equations, 



In conclusion, we derived exact hydrodynamic equations 
from the linearized Boltzmann-BGK equation. The main nov- 
elty is the numerical non-perturbative procedure to solve the 
invariance equation. In turn, the highly efficient numerical 
approach is made possible by choosing a convenient coor- 
dinate system and establishing symmetries of the invariance 
equation. The invariant manifold in the space of distribution 
functions is thereby completely characterized, that is, not only 
equations of hydrodynamics are obtained but also the corre- 
sponding distribution function is made available. The perti- 
nent data can be used, in particular, as a much needed bench- 
mark for computation-oriented kinetic theories such as lattice 
Boltzmann models, as well as for constructing novel models 
using quadratures in the velocity space [13, 14]. Finally, we 
have established a novel non-perturbative computational ap- 
proach to finding invariant manifolds of kinetic equations. 

The present approach can be extended to the Boltzmann 
equation with other collision terms. The above derivation of 
hydrodynamics is done under the standard assumption of lo- 
cal equilibrium [1], however the assumption itself is open to 
further study [15] [We thank H.C. Ottinger for this important 
remark]. I.V.K. acknowledges support of CCEM-CH. 
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